function idat = fsb_remove_outliers(idat,imean)

% FSB - DEV: remove outliers from image data
%
% EXAMPLE:
% idat = fsb_remove_outliers(idat,imean)
%
% INPUT:
% idat:     4D image data
% imean:    Mean image data
%
% OUTPUT:
% idat:     Corrected 4D image data
% 
% CALLED BY:
% FSB.m
%
% NOTES:
% In development.
% Outliers are defined as datapoints in a
% voxel time course that are more than 3
% standard deviations above or below the mean
% voxel value over time.
% Outliers are then set to either 3 SD above
% or below the mean voxel time course value.
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 0.1
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Outlier correction - already optimized for speed
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Removing outliers...');
idat = single(idat);
idat(idat==0)=NaN;
idat(idat>32000) = NaN;
waitbar(10/100);
mean_idat = single(imean);
waitbar(20/100);
mean_idat = repmat(mean_idat,[1 1 1 size(idat,4)]);
waitbar(50/100);
std_idat = fsb_calculate_std(idat);
std_idat2 = repmat(std_idat,[1 1 1 size(idat,4)]);
waitbar(60/100);
idat_log1 = single(idat)>(mean_idat+(3*std_idat2));
waitbar(80/100);
idat_log2 = single(idat)<(mean_idat-(3*std_idat2));
waitbar(90/100);
max_idat = int16(mean_idat+(3*std_idat2));
waitbar(95/100);
min_idat = int16(mean_idat-(3*std_idat2));

if min_idat<1
    min_idat = 0;
end

waitbar(100/100);
idat(idat_log1)= max_idat(idat_log1);
idat(idat_log2)= min_idat(idat_log2);
disp(['Upper outlier correction threshold: ' num2str(max(max_idat(:)))]);
disp(['Lower outlier correction threshold: ' num2str(min(min_idat(:)))]);
idat(idat<0)=0;
idat = int16(idat);
close(h);

end
